#include <MapKit/MapKit.h>
#include <cstdlib>
#include <queue>
#include <glm/glm.hpp>
#include <glm/gtc/matrix_transform.hpp>

void MapKit::Randomize(const unsigned int& seed)
{
    srand(seed);
}

double MapKit::RandDouble()
{
    union
    {
        double_t d;
        uint64_t u64;
        uint16_t u16[sizeof(uint64_t) / sizeof(uint16_t)];
    } x;
    for (uint32_t i = 0; i < sizeof(uint64_t) / sizeof(uint16_t); ++i)
        x.u16[i] = (rand() << 1) | (rand() & 1);
    x.u64 = (x.u64 >> 12) | 0x3FF0000000000000ULL;
    return x.d - 1;
}

void MapKit::GetPoisson(const double& r, const uint32_t& k)
{
    // according to https://blog.csdn.net/gouki04/article/details/100143747
    // step 0
    uint32_t n = 2;
    double cellSize = r / glm::sqrt(n);
    uint32_t cols = (uint32_t)glm::ceil(width / cellSize);
    uint32_t rows = (uint32_t)glm::ceil(height / cellSize);
    std::vector<glm::dvec2> cells;
    std::vector<std::vector<int>> grids(rows, std::vector<int>(cols, -1));
    // step 1
    auto v0 = glm::dvec2(RandDouble() * width, RandDouble() * height);
    uint32_t col = (uint32_t)glm::floor(v0.x / cellSize);
    uint32_t row = (uint32_t)glm::floor(v0.y / cellSize);
    int32_t v0Idx = cells.size();
    cells.push_back(v0);
    grids[row][col] = (int32_t)v0Idx;
    std::vector<int32_t> activeList;
    activeList.push_back(v0Idx);
    // step 2
    while (activeList.size() > 0)
    {
        auto rid = rand() % activeList.size();
        auto viIdx = activeList[rid];
        auto& vi = cells[viIdx];
        bool found = false;
        for (uint32_t i = 0; i < k; i++)
        {
        TryFindNewVertex:
            auto dir = [&]()->glm::dvec2 {
                glm::dvec2 rd;
                double r1 = r, r2 = 2 * r;
                double a = 1 / (r2 * r2 - r1 * r1);
                double rr = glm::sqrt(RandDouble() / a + r1 * r1);
                double rt = RandDouble() * glm::two_pi<double>();
                rd.x = rr * glm::cos(rt);
                rd.y = rr * glm::sin(rt);
                return rd;
            }();
            auto vk = vi + dir;
            if (vk.x < 0 || vk.x >= width || vk.y < 0 || vk.y >= height)
                goto TryFindNewVertex;
            col = (uint32_t)glm::floor(vk.x / cellSize);
            row = (uint32_t)glm::floor(vk.y / cellSize);
            if (grids[row][col] != -1)
                continue;
            bool ok = true;
            int32_t minR = (int32_t)glm::floor((vk.y - r) / cellSize);
            int32_t maxR = (int32_t)glm::floor((vk.y + r) / cellSize);
            int32_t minC = (int32_t)glm::floor((vk.x - r) / cellSize);
            int32_t maxC = (int32_t)glm::floor((vk.x + r) / cellSize);
            for (int32_t oR = minR; oR <= maxR; oR++)
            {
                if (oR < 0 || oR >= (int32_t)rows)
                    continue;
                for (int32_t oC = minC; oC <= maxC; oC++)
                {
                    if (oC < 0 || oC >= (int32_t)cols)
                        continue;
                    auto vjIdx = grids[oR][oC];
                    if (vjIdx != -1)
                    {
                        auto& vj = cells[vjIdx];
                        auto dist = glm::length(vj - vk);
                        if (dist < r)
                        {
                            ok = false;
                            goto EndOfDistanceCheck;
                        }
                    }
                }
            }
        EndOfDistanceCheck:
            if (ok)
            {
                auto vkIdx = cells.size();
                cells.push_back(vk);
                grids[row][col] = vkIdx;
                activeList.push_back(vkIdx);
                found = true;
                break;
            }
        }
        if (!found)
            for (auto it = activeList.begin(); it != activeList.end();)
                if ((*it) == viIdx)
                {
                    it = activeList.erase(it);
                    break;
                }
                else
                    it++;
    }
    for (auto p : cells)
    {
        coords->push_back({ p.x, p.y });
        MKVertex* vt = MKVertex::Create();
        vt->x = p.x;
        vt->y = p.y;
        bool found = false;
        for (MKVertex* q : *vertices)
            if (*q == *vt)
            {
                found = true;
                break;
            }
        if (!found)
            vertices->push_back(vt);
    }
}

void MapKit::FractalAltitude(double ratio)
{
    int rw = (int)(glm::ceil(glm::log(width) / glm::log(2)));
    int rh = (int)(glm::ceil(glm::log(height) / glm::log(2)));
    int ft = std::max(rw, rh);
    int ff = (1 << ft) + 1;
    fractalSize = ff;
    std::vector<std::vector<double>> board(ff, std::vector<double>(ff, 0));
    for (int t = 0; t < ft; t++)
    {
        // diamond
        int sxy = ff >> (t + 1);
        int count = 1 << t;
        int step = 1 << (ft - t);
        double rough = glm::pow(ratio, t);
        for (int i = 0; i < count; i++)
            for (int j = 0; j < count; j++)
            {
                double r = (RandDouble() * 2 - 1) * rough;
                int cx = sxy + i * step, cy = sxy + j * step;
                board[cx][cy] = (
                    board[cx - (1 << (ft - t - 1))][cy - (1 << (ft - t - 1))] +
                    board[cx + (1 << (ft - t - 1))][cy - (1 << (ft - t - 1))] +
                    board[cx - (1 << (ft - t - 1))][cy + (1 << (ft - t - 1))] +
                    board[cx + (1 << (ft - t - 1))][cy + (1 << (ft - t - 1))]
                ) / 4 + r;
            }
        // square
        int quad = 1 << (ft - t - 1);
        int amount = (1 << (t + 1)) + 1;
        for (int i = 0; i < amount; i++)
            for (int j = 0; j < amount; j++)
            {
                if ((i + j) % 2 == 1)
                {
                    double r = (RandDouble() * 2 - 1) * rough;
                    int cx = i * quad, cy = j * quad;
                    board[cx][cy] = (
                        board[i == 0 ? ff - 1 - quad : cx - quad][cy] +
                        board[cx][j == 0 ? ff - 1 - quad : cy - quad] +
                        board[i * quad == ff - 1 ? quad : cx + quad][cy] +
                        board[cx][j * quad == ff - 1 ? quad : cy + quad]
                    ) / 4 + r;
                }
            }
    }
    double value = glm::abs(board[0][0]);
    for (int i = 0; i < ff; i++)
        for (int j = 0; j < ff; j++)
            if (glm::abs(board[i][j]) > value)
                value = glm::abs(board[i][j]);
    for (int i = 0; i < ff; i++)
        for (int j = 0; j < ff; j++)
            board[i][j] /= value;
    *altitude = board;
}

void MapKit::PerlinNoiseAltitude()
{
    int rw = (int)(glm::ceil(glm::log(width) / glm::log(2)));
    int rh = (int)(glm::ceil(glm::log(height) / glm::log(2)));
    int ft = std::max(rw, rh);
    int ff = (1 << ft) + 1;
    fractalSize = ff;
    std::vector<std::vector<double>> board(ff, std::vector<double>(ff, 0));
    *altitude = board;
}

double MapKit::AltitudeLinearSampler2D(const double& x, const double& y)
{
    int x0 = int(glm::floor(x)), y0 = int(glm::floor(y));
    int x1 = int(glm::ceil(x)), y1 = int(glm::ceil(y));
    double d = (x1 - x0) * (y1 - y0);
    if (d == 0)
        return (*altitude)[x0][y0];
    double q00 = (*altitude)[x0][y0];
    double q10 = (*altitude)[x1][y0];
    double q01 = (*altitude)[x0][y1];
    double q11 = (*altitude)[x1][y1];
    return q00 / d * (x1 - x) * (y1 - y) + q10 / d * (x - x0) * (y1 - y) + q01 / d * (x1 - x) * (y - y0) + q11 / d * (x - x0) * (y - y0);
}